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1 Introduction 

Many solids, such as alloys or doped semiconductors, form crystals consisting 
of two or more chemical species. In contrast to amorphous systems they show 
a regular, periodic arrangement of sites but the different species are statisti- 
cally distributed over the available lattice sites. This type of disorder is often 
called compositional. Likewise crystal imperfections present in any real mate- 
rial give rise to substitutional disordered structures. The presence of disorder 
has drastic effects on various characteristics of physical systems, most notably 
on their electronic structure and transport properties. 

A particularly interesting case is a lattice composed of accessible and in- 
accessible sites (bonds). Dilution of the accessible sites defines a percolation 
problem for the lattice which can undergo a geometric phase transition be- 
tween a connected and a disconnected phase. Since 'absent' sites (bonds) act 
as infinite barriers such a model can be used to describe percolative particle 
transport through random resistor networks (see Fig. 1). Another example is 
the destruction of magnetic order in diluted classical magnets. The central 
question of classical percolation theory is whether the diluted lattice contains 
a cluster of accessible sites that spans the entire sample in the thermodynamic 
limit or whether it decomposes into small disconnected pieces. 

The corresponding problem of percolation of a Quantum particle through 
a random medium is much more involved. Here the famous concept of An- 
derson localization comes into play [1] . Anderson argued that the one-particle 
wave functions in macroscopic, disordered quantum systems at zero tempera- 
ture can be exponentially localized. This happens if the disorder is sufficiently 
strong or in energy regions where the density of states is sufficiently small [2] . 
The transition from extended to localized states is a result of quantum inter- 
ference arising from elastic scattering of the particle at the randomly varying 
site energies. Particles that occupy exponentially localized states are restricted 
to finite regions of space. On the other hand, particles in extended states can 
escape to infinity and therefore contribute to transport (see Fig. 2). 
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Fig. 1. 2D site percolation. Shown are lattice realizations below (left) and above 
(right) the classical percolation threshold p c — 0.592746 



For the quantum percolation problem all one-particle states are localized 
below the classical percolation threshold p c , of course. Furthermore for arbi- 
trary concentrations of conducting sites p > p c there always exist localized 
states on the infinite spanning cluster [3] . Since for a completely ordered sys- 
tem (p = 1) all states are extended and no states are extended for p < p c , one 
might expect a disorder-induced Anderson transition at some critical concen- 
tration p q , with p c <p q < 1. The reason is that the quantum nature of parti- 
cles makes it harder for them to pass through narrow channels, despite the fact 
that quantum particles may tunnel through classically forbidden regions [4]. 
As yet the existence of a quantum percolation threshold, different from the 
classical one, is discussed controversial, in particular for the two-dimensional 
(2D) case where also weak localization effects might become important [5]. 

Localization and quantum percolation have been the subject of much re- 
cent research; for a review of the more back dated work see e.g. [6, 7). The 
underlying Anderson model (with uniformly distributed local potentials) and 
site percolation or binary alloy models (with a bimodal distribution of the 
on-site energies) are the two standard models for disordered solids. Although 
both problems have much in common, they differ in the type of disorder, and 
the localization phenomena of electrons in substitutional alloys are found to 
be much richer than previously claimed. For instance, the binary alloy model 
exhibits not one, as the 31? Anderson model, but several pairs of mobility 
edges, separating localized from extended states [3, 8, 9]. Moreover it appears 
that 'forbidden energies' exist, in the sense that near the band center the 
density of states continuously goes to zero. These effects might be observed 
in actual experiments. 

Understanding these issues is an important task which we will address in 
this tutorial-style review by a stochastic approach that allows for a compre- 
hensive description of disorder on a microscopic scale. This approach, which 
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Fig. 2. Anderson transition in disordered systems 

we call the local distribution (LD) approach [9, 10], considers the local density 
of states (LDOS), which is a quantity of primary importance in systems with 
prominent local interactions or scattering. What makes the LD approach 'non- 
standard' is that it directly deals with the distribution of the LDOS in the 
spirit that Anderson introduced in his pioneering work [1]. While the LDOS 
is just related to the amplitude of the electron's wave function on a certain 
lattice site, its distribution captures the fluctuations of the wave function 
through the system. 

The present paper will be organized as follows. After introducing the basic 
concepts of the LD approach in Sect. 2 and explaining how to calculate the 
LDOS via the highly efficient Kernel Polynomial Method (KPM), we exem- 
plify the LD concept by a study of localization within the Anderson model. 
Then, having all necessary tools at hand, we proceed in Sect. 3 to the prob- 
lem of quantum percolation. Refining previous studies of localization effects in 
quantum percolation, we will demonstrate the 'fragmentation' of the spectrum 
into extended and localized states. The existence of a quantum percolation 
threshold is critically examined. To this end, we investigate for the 2D site- 
percolation problem the dynamical properties of an initially localized wave 
packet on the spanning cluster, using Chebyshev expansion in the time do- 
main. In Sect. 4 we apply the findings to several classes of advanced materials 
where transport is related to percolating current patterns. Here prominent 
examples are 2D undoped graphene and the 3D colossal magnetoresistive 
manganites. Finally, in Sect. 5 we conclude with a short summary of the 
topic. 

2 Local Distribution Approach 
2.1 Conceptional background 

In the theoretical investigation of disordered systems it turned out that the 
distribution of random quantities takes the center stage [1, 11]. Starting from 
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a standard tight-binding Hamiltonian of noninteracting electrons on a D- 
dimensional hypcr-cubic lattice with TV = L D sites, the effect of compositional 
disorder in a solid may be described by site-dependent local potentials ej . 

N 

H = -* E( c h- + H - c -) + E e A c 3 ■ (!) 

The operators cj (cj) create (annihilate) an electron in a Wannier state cen- 
tered at site j, the on-site potentials cj are drawn from some probability 
distribution p[ej], and t denotes the nearest-neighbor hopping element. While 
all characteristics of a certain material are determined by the corresponding 
distribution p[ej], we have to keep in mind, that each actual sample of this 
material constitutes only one particular realization {ej}. The disorder in the 
potential landscape breaks translational invariancc, which normally can be 
exploited for the description of ordered systems. Hence, we have to focus on 
sitc-dependent quantities like the LDOS at lattice site i, 

N 

Pi (E) = J2\(i\n)\ 2 S(E-E n ), (2) 

71=1 

where \n) is an eigenstate of H with energy E n and \i) — c||0). Probing differ- 
ent sites in the crystal, Pi{E) will strongly fluctuate, and recording the values 
of the LDOS we may accumulate the distribution f[pi(E)]. In the thermo- 
dynamic limit, taking into account infinitely many lattice sites, the shape of 
f[pi(E)] will be independent of the actual realization of the on-site potentials 
{e^} and the chosen sites, but depend solely on the underlying distribution 
p[ej]. Thus, in the sense of distributions, we have restored translational invari- 
ance, and the study of f[pi{E)\ allows us to discuss the complete properties 
of the Hamiltonian (1). The probability density f[pi(E)] was found to have 
essentially different properties for extended and localized states [12]. For an 
extended state the amplitude of the wave function is more or less uniform. Ac- 
cordingly f[pi(E)} is sharply peaked and symmetric about the (arithmetically 
averaged) mean value p mc (E). On the other hand, if states become localized, 
the wave function has considerable weight only on a few sites. In this case 
the LDOS strongly fluctuates throughout the lattice and the corresponding 
LDOS distribution is mainly concentrated at pi — 0, very asymmetric and has 
a long tail. The rare but large LDOS-values dominate the mean value p mc (E), 
which therefore cannot be taken as a good approximation of the most prob- 
able value of the LDOS. Such systems are referred to as 'non-self-averaging'. 
In numerical calculations, this different behavior may be used to discriminate 
localized from extended states in the following manner: Using the KPM with 
a variable number of moments for different interval sections (see Sect. 2.2) we 
calculate the LDOS for a large number of samples, K r , and sites, K s . Note 
that from a conceptual point of view, it makes no difference if we calculate 
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f[pi(E)] for one large sample and many sites or consider smaller (but not too 
small) systems and more realizations of disorder. As the latter procedure is 
numerically more favorable, we will revert to this. Next, we sort the LDOS 
values for all energies within a small window around the desired energy E into 
a histogram. As the LDOS values vary over several orders of magnitude, a log- 
arithmic partition of the histogram presents itself. Although this procedure 
allows for the most intuitive determination of the localization properties in the 
sense of Anderson [1], there are two drawbacks. First, to achieve a reasonable 
statistics, in particular for the slots at small LDOS values, a huge number 
of realizations is necessary. To alleviate the problems of statistical noise, it is 
advantageous to look at integral quantities like the distribution function 

F[p i (E)}= I*™ ftfWWAE) , (3) 
Jo 

which also allows for the determination of the localization properties. While 
for extended states the more or less uniform amplitudes lead to a steep rise 
of F[pi], for localized states the increase extends over several orders of magni- 
tude. Second, for practical calculations the recording of the whole distribution 
(or distribution function) is a bit inconvenient, especially if we want to dis- 
cuss the localization properties of the whole band. Therefore we try to capture 
the main features of the distribution by comparing two of its characteristic 
quantities, the (arithmetically averaged) mean DOS, 

r s fe=li=l 

and the (geometrically averaged) so called 'typical' DOS , 

Pty (E) = cxp — — J2 E HP*( E )) ■ (5) 
^ r s fc=ii=i ' 

The typical DOS puts sufficient weight on small values of p,. Therefore com- 
paring pty(E) and p me (E) allows to detect the localization transition. We 
classify a state at energy E with p me (E) ^ as localized if pt y (E) = and 
as extended if Pt y {E) ^ 0. This method has been applied successfully to the 
pure Anderson model [13, 14, 15] and even to more complex situations, where 
the effects of correlated disorder [16], electron-electron interaction [17, 18] or 
electron-phonon coupling [19, 20] were taken into account. 

2.2 Calculation of the Local Density of States by the 
Variable Moment Kernel Polynomial Method 

At first glance, (2) suggests that the calculation of the LDOS could require a 
complete diagonalization of H . It turns out, however, that an expansion of pi 
in terms of a finite series of Chebyshev polynomials T n (x) = cos(n arccosrr) 
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allows for an incredibly precise approximation [21]. Since the Chebyshev poly- 
nomials form an orthogonal set on [—1,1], prior to an expansion the Hamilto- 
nian H needs to be rescaled, H = (H — b)/a. For reasons of numerical stability, 
we choose the parameters a and b such that the extreme eigenvalues of H are 
±0.99. In this way, the outer parts of the interval, where the strong oscilla- 
tions of T n (x) can amplify numerical errors, contain no physical information 
and may be discarded. In terms of the coefficients, the so called Chebyshev 
moments, 

} 

Mm = J Pi(x)T m {x)dx = ^2(i\n)(n\T m (x n )\i) = (i\T m (H)\i) , (6) 



the approximate LDOS reads 



M-l 



Pi{x) = x 2 {doVo + 2 9m^ m T m {x) J • (7) 



m.—l 



The Gibbs damping factors 

(TYl(f) \ (f) 
1 cos(m</>) H — sin(m</<) cot(0) (8) 
7T / 7T 

with (f> — ir/(M + 1), are introduced to damp out the Gibbs oscillations 
resulting from finite-order polynomial approximations. The introduction of 
these factors corresponds to convoluting the finite series with the so called 
Jackson kernel [22]. In essence, a 5-peak at position xo is approximated by an 
almost Gaussian of width [211 



g = V 2(M + D ^-^-Mf-^ir- ^ 

Thus for a fixed number of moments, M, the resolution of the expansion 
gets better towards the interval boundaries. While in most applications this 
feature is harmless or even useful, here a uniform resolution throughout the 
whole band is mandatory to discriminate resolution effects from localization. 
This gets clear, if we respect that a small value of pi at a certain energy 
may either be due to a true low amplitude of the wave function, or to the 
absence of any energy level for the current disorder realization within one 
kernel width. Depending on which part of the interval we want to reconstruct, 
we need to restrict the used number of moments in (7) accordingly to ensure 
a constant cr. We call this procedure the Variable Moment KPM (VMKPM). 
The resulting approximations of a series of 5-peaks using the standard KPM 
and the VMKPM, respectively, are compared in Fig. 3. 

For practical calculation of the moments, we may profit from the recursion 
relations of the Chebyshev polynomials, 
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Fig. 3. Approximation of a series of equally weighted 5-peaks using the standard 
KPM (left) and the VMKPM (right). While both methods reproduce the correct 
spectral weight (solid line), only in the latter case a uniform resolution is obtained, 
reflected by the equal height of the peaks. The dashed-dotted line in the left panel 
is a guide to the eye being proportional to the inverse of the resolution (9) 



T m+1 {x) = 2xT m (x) - T m -i[x) , (10) 

starting from Tq(x) = 1 and T\(x) = x, and calculate the \i m iteratively. 
Additionally, we may reduce the numerical effort by another factor 1/2 by 
generating two moments with each matrix vector multiplication by H, 

A*2m-i = 2(i\T m (H)T m -i(H)\i) - n\ , 
[L2m = 2(i\T m (H)T m (H)\i) - ^ • 

Note that the algorithm requires storage only for the sparse matrix H and two 
vectors of the corresponding dimension. Having calculated the desired num- 
ber of moments, we calculate several reconstructions (7) for different M. We 
obtain the final result with uniform resolution by smoothly joining the corre- 
sponding results for the different subintervals. As the calculation of the [i m 
dominates the computational effort, the additional overhead for performing 
several reconstructions is negligible as they can be done using a Fast Fourier 
Transform (FFT) routine. 



2.3 Illustration of the Method: Anderson Localization in 3D 

To clarify the power of the method, let us briefly apply it to the Anderson 
model of localization [1], for which the principal results are well known and 
can be found in the literature [2, 7]. The Anderson model is described by the 
Hamiltonian (1), using local potentials ej, which are assumed to be indepen- 
dent, uniformly distributed random variables, 

P[^] = ^Q-\e j \) , (12) 

where the parameter 7 is a measure for the strength of disorder. The spectral 
properties of the Anderson model have been carefully analyzed (see e.g. [23]). 
For sufficiently large disorder or near the band tails, the spectrum consists 
exclusively of discrete eigenvalues, and the corresponding eigenfunctions are 
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Fig. 4. Upper panel: Finite-size scaling of the normalized LDOS distribution 
/[Pi/Pmo] at the band center (\E\ < O.Olt) for the Anderson model. For the dif- 
ferent system sizes TV we adapted the resolution a such that TVer ~ 6.14 is con- 
stant and calculated 10 6 , 10 4 , 1024, 400 realizations of p t for TV = 10 3 , 20 3 , 40 3 , 80 3 . 
Lower left panel: Double- logarithmic plot of /[p;/pme] for the localized case to- 
gether with log-normal fits (13) to the data. Note that for better visibility the data 
for TV = 40 3 ,20 3 , 10 3 were shifted vertically by 2, 3, 4 orders of magnitude towards 
smaller values. Right lower panel: Distribution function Ffpi/pme] of the above data 



exponentially localized. Since localized electrons do not contribute to the 
transport of charge or energy, the energy that separates localized and ex- 
tended eigenstates is called the mobility edge. For any finite amount of dis- 
order 7 > 0, on a II? lattice, all eigenstates of the Anderson Hamiltonian 
are localized [24, 25]. This is believed to hold also in 2D, where the existence 
of a transition from localized to extended states at finite 7 would contradict 
the one parameter scaling theory [26, 27]. In three dimensions, the disorder 
strength has a more distinctive effect on the spectrum. Only above a critical 
value 7 C all states are localized, whereas for 7 < 7 C a pair of mobility edges 
separates the extended states in the band center from the localized ones near 
the band edges [28]. For this reason, the 3D case serves as a prime example 
on which we demonstrate how to discriminate localized from extended states 
within the local distribution approach. In the upper panel of Fig. 4 we show 
the resulting distribution of Pi{E — 0), normalized by its mean value p mc , 
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Fig. 5. Mean (solid) and typical (dashed) density of states of the 3D Anderson 
model on a lattice with N = 60 3 sites and periodic boundary conditions (PBC) 
for different values of disorder 7. We used Na = 45 and the shown curves result 
from an ensemble of 2000 individual LDOS spectra obtained for different sites and 
realizations of disorder. Lower right panel: Mobility edge (dots) as obtained for 
(Pty/pme)c = 0.05 and Lifshitz boundaries (dashed line) 



for two characteristic values of disorder. As /9 me is a function of disorder, this 
normalization ensures (pi/p mc ) = 1 independent of 7, allowing for an appro- 
priate comparison. In the delocalized phase, 7 — 3t, the distribution is rather 
symmetric and peaked close to its mean value. Note that increasing both the 
system size and VMKPM resolution, such that the ratio of mean level spacing 
to the width of the Jackson kernel is fixed, does not change the distribution. 
This is in strong contrast to the localized phase, 7 = 24t, where the distribu- 
tion of Pi(E) is extremely asymmetric. Although most of the weight is now 
concentrated close to zero, the distribution extends to very large values of 
Pi, causing the mean value to be much larger than the most probable value. 
Performing a similar finite-size scaling underlines both the asymmetry and 
the singular behavior which we expect for infinite resolution in the thermo- 
dynamic limit. Note also, that for the localized case, the distribution of the 
LDOS is extremely well approximated by a log- normal distribution [29], 

. . 11 / {\n(xlx )f\ 

^-^i^l-Hsr-J' (13) 

as illustrated in the lower left panel of Fig. 4. The shifting of the distribution 
towards zero for localized states is most obvious in the distribution function 
F[pi] which is depicted in the lower right panel of Fig. 4. While for extended 
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states the more or less uniform amplitudes lead to a steep rise of F[pi], for 
localized states the increase extends over several orders of magnitude. Cap- 
turing the essential features of the LDOS distribution by concentrating on the 
mean (p mc ) and typical (p ty ) density of states, we determine the localization 
properties for the whole energy band depending on the disorder. As can be 
seen from Fig. 5, p mc {E) and pty(E) are almost equal for extended states, 
whereas for localized states pty{E) vanishes while p mc (E) remains finite. Us- 
ing the well-established value ^ciE = 0) ~ 16. 5t as a calibration for the 
critical ratio pty/pmc, required to distinguish localized from extended states 
for the used system size and resolution, we reproduce the mobility edge in the 
energy-disorder plane [30, 2, 31] (see lower right panel of Fig. 5). We also find 
the well-known re-entrant behavior near the unperturbed band edges [32, 33]: 
Varying 7 for some fixed values of E (6t < E < 7.6i) a region of extended 
states separates two regions of localized states. The Lifshitz boundaries, shown 
as dashed lines, indicate the energy range, where eigenstates are in principle 
allowed. As the probability of reaching the Lifshitz boundaries is exponen- 
tially small, we cannot expect to find states near these boundaries for the 
finite ensembles considered in any numerical calculation. 

3 Localization Effects in Quantum Percolation 

In disordered solids the percolation problem is characterized by the interplay 
of pure classical and quantum effects. Besides the question of finding a perco- 
lating path of 'accessible' sites through a given lattice, the quantum nature of 
the electrons imposes further restrictions on the existence of extended states 
and, consequently, of a finite DC-conductivity. 

As a particularly simple model we start again from the basic Hamilto- 
nian (1) drawing the {ej} from the bimodal distribution 

Pfo] =p8(e j -e A ) + {l-p)5(e j -e B ) . (14) 

The two energies ca and cb could, for instance, represent the potential land- 
scape of a binary alloy A p Bi_ p , where each site is occupied by an A or B atom 
with probability por 1 — p, respectively. Therefore we call (1) together with 
the distribution (14) the binary alloy model. In the limit A = (cb — za) — * 00 
the wave function of the A sub-band vanishes identically on the -B-sites, mak- 
ing them completely inaccessible for the quantum particles. We then arrive at 
a situation where non-interacting electrons move on a random ensemble of iV 
lattice points, which, depending on p, may span the entire lattice or not. The 
corresponding quantum site percolation model reads 

H = -tJ2 (ctcj+H.c.) , (15) 

(ij)eA 

where the summation extends over nearest-neighbor A-sites only and, without 
loss of generality, €a is chosen to be zero. 
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Within the classical percolation scenario the percolation threshold p c is 
defined by the occurrence of an infinite cluster of connected A sites. For 
the simple cubic lattice this site-percolation threshold is p^ c D = 0.311609 [34] 
for the 31? case and p^P = 0.592746 [35] in 2D. In the quantum case, the 
multiple scattering of the particles at the irregular boundaries of the cluster 
can suppress the wave function, in particular within narrow channels or close 
to dead ends of the cluster. Hence, this type of disorder can lead to absence 
of diffusion due to localization, even if there is a classical percolating path 
through the crystal. On the other hand, for finite A the tunneling between A 
and B sites may cause a finite DC-conductivity although the A sites are not 
percolating. Naturally, the question arises whether the quantum percolation 
threshold p q , given by the probability above which an extended wave function 
exists within the A sub-band, is larger or smaller than p c . Previous results [36, 
37] for finite values of A indicate that the tunneling effect has a marginal 
influence on the percolation threshold as soon as A >• AtD. 

3.1 3D Site Percolation 

Before discussing possible localization phenomena let us investigate the be- 
havior of the mean DOS for the binary alloy and quantum percolation model 
in 3D. Figure 6 shows that as long as €a and es do not differ too much 
there exists a single and if p ^ 0.5 asymmetric electronic band [36]. At about 
A ~ AtD this band separates into two sub-bands centered at ca and e^, re- 
spectively. The most prominent feature in the split-band regime is the series 
of spikes at discrete energies within the band. As an obvious guess, we might 
attribute these spikes to eigenstates on islands of A or B sites being isolated 
from the main cluster [3, 38]. It turns out, however, that some of the spikes 
persist, even if we neglect all finite clusters and restrict the calculation to 
the N sites of A-type on the spanning cluster, Aoq. This is illustrated in the 
upper panels of Fig. 7, where we compare the DOS of the binary alloy model 
at A — > oo and the quantum percolation model. Increasing the concentration 
of accessible sites the mean DOS of the spanning cluster is evocative of the 
DOS of the simple cubic lattice, but even at large values of p a sharp peak 
structure remains at E = (see Fig. 7, lower panels). 

To elucidate this effect, which for a long time was partially not accounted 
for in the literature [3, 36, 39, 40], in more detail, in Fig. 8 we fix p at 
0.337, shortly above the classical percolation threshold. In addition to the 
most dominant peaks at E/t = 0, ±l,±\/2, we can resolve distinct spikes 
at E/t = \ (±1 ± \/5) , ±V3,±V2±72, . . . . These special energies coincide 
with the eigenvalues of the tight-binding model on small clusters of the ge- 
ometries shown in the right part of Fig. 8. In accordance with [3] and [41] 
we can thus argue that the wave functions, which correspond to these special 
energies, are localized on some 'dead ends' of the spanning cluster. 

The assumption that the distinct peaks correspond to localized wave func- 
tions is corroborated by the fact that the typical DOS vanishes or, at least, 
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Fig. 6. Mean (upper solid line) and typical (lower dashed line) DOS of the binary 
alloy model on a N = 64 3 lattice with PBC for different concentrations of ^4-sites p 
and local potential differences A. Taking all sites into account, we chose the resolu- 
tion such that Na = 45 and calculated 1000 individual LDOS spectra for different 
probe realizations and sites 

shows a dip at these energies. Occurring also for finite A (Fig. 6), this effect 
becomes more pronounced as A — ► oo and in the vicinity of the classical per- 
colation threshold p c . From the study of the Anderson model we know that 
localization leads at first to a narrowing of the energy window containing ex- 
tended states. For the percolation problem, in contrast, with decreasing p the 
typical DOS indicates both localization from the band edges and localization 
at particular energies within the band. Since finite cluster wave functions like 
those shown in Fig. 8 can be constructed for numerous other, less probable 
geometries [42], an infinite discrete series of such spikes might exist within 
the spectrum, as proposed in [41]. The picture of localization in the quantum 
percolation model is then quite remarkable. If we generalize our numerical 
data for the peaks at E — and E — ±t, it seems as if there is an infinite 
discrete set of energies with localized wave functions, which is dense within 
the entire spectrum. In between there are many extended states, but to avoid 
mixing, their density goes to zero close to the localized states. Facilitated by 
its large special weight (up to 11% close to p c ) this is clearly observable for 
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Fig. 7. Mean (solid line) and typical (dashed line) DOS for the quantum percolation 
model in the limit A — ► oo. While in the upper left panel all A-sites are taken into 
account, the other three panels show data for the restricted model on the spanning 
cluster Aqo only (note that p ty is smaller in the former case because there are more 
sites with vanishing amplitude of the wave function). Using a — 0.0002, we adapt 
the system sizes to ensure Na = 45 for the A^, case, i.e. , N = 64 3 ,70 3 ,84 3 for 
p = 0.86, 0.66, 0.413. For the case of all A sites at p = 0.413, N = 82 3 guarantees 
Na = 45. For each system, we consider 512 individual LDOS spectra for different 
sites and realizations of disorder 
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Fig. 8. Left: Mean (solid line) and typical (dashed line) DOS for the quantum 
percolation model with p — 0.337 on a TV = 100 3 lattice (restriction to A^,, PBC, 
Na = 45). The data are an ensemble average over 1000 realizations of disorder and 
different sites. Note that p ty < 10~ 3 holds in the whole band. Right: Some cluster 
configurations related to the special energies at which the peaks in p mo occur 
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the peak at E = 0, and we suspect similar behavior at E = ±t. For the other 
discrete spikes the resolution of our numerical data is still too poor and the 
system size might be even too small to draw a definite conclusion. 





Fig. 9. Normalized probability density 7V|(i|n)| 2 for three eigenstates of the quan- 
tum percolation model. The chosen energies are representative for a localized 
(E — 3.501i), extended (E — 0.64SM) and anomalous localized (E = 0) state. 
We performed an exact diagonalization for Aoc for one realization of disorder on 
a N = 32 3 lattice with PBC and p = 0.458 



In order to understand the internal structure of the extended and localized 
states we calculate the probability density of specific eigenstates of (15) re- 
stricted to Aoo for a fixed realization of disorder. Figure 9 visualizes the spatial 
variation of |(i|n)| 2 for an occupation probability well above the classical per- 
colation threshold. The figure clearly indicates that the state with E = 0.649i 
is extended, i.e. the spanning cluster is quantum mechanically 'transparent'. 
On the contrary, at E = 3.501t, the wave function is completely localized on 
a finite region of the spanning cluster. This is caused by the scattering of the 
particle at the random surface of the spanning cluster. 

A particularly interesting behavior is observed at E = 0. Here, the eigen- 
states are highly degenerate and we can form wave functions that span the 
entire lattice in a checkerboard structure with zero and non-zero amplitudes 
(see Fig. 9). Although these states are extended in the sense that they are not 
confined to some region of the cluster, they are localized in the sense that they 
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do not contribute to the DC-conductivity. This is caused by the alternating 
structure which suppresses the nearest-neighbor hopping, and in spite of the 
high degeneracy, the current matrix element between different E = states is 
zero. Hence, having properties of both classes of states these states are called 
anomalously localized [43, 44]. Another indication for the robustness of this 
feature is its persistence for mismatching boundary conditions, e.g., periodic 
(anti-periodic) boundary conditions for odd (even) values of the linear exten- 
sion L. In these cases the checkerboard is matched to itself by a manifold of 
sites with vanishing amplitude. 

Furthermore for the previously mentioned special energies, the wave func- 
tion vanishes identically except for some finite domains on loose ends (like 
those shown in the right panel of Fig. 8), where it takes, except for normali- 
sation, the values (±1, ^=1), (1, ±\/2, 1), (— 1, 0, 1), ... Note that these regions 
are part of the spanning cluster, connected to the rest of Aoo by sites with 
wave function amplitude zero [41] . 




Fig. 10. Left panel: Normalized typical DOS pty/pmc in the concentration-energy 
plane for the quantum percolation model restricted to A^. For a fixed resolu- 
tion of a = 0.0002 we adapted the system sizes to keep Na — 45 constant, i.e., 
N = 62 3 , . . . , 120 3 for p = 0.944, . . . , 0.318. The results base on the calculation 
of 1000 individual LDOS spectra for different cluster realizations and sites. Right 
panel: Distribution function of the LDOS for four characteristic states, for which the 
parameters (p, E) are indicated by crosses in the left panel. The parameters for the 
curves A, B, C are the same for which we showed the characteristic states in Fig. 9 

In the past most of the methods used in numerical studies of Anderson 
localization have also been applied to the binary alloy model and the quan- 
tum percolation model in order to determine the quantum percolation thresh- 
old pg, defined as the probability p below which all states are localized (see, 
e.g., [45, 46] and references therein). The existence of p q is still disputed. As 
yet the results for p q are far less precise than, e.g., the values of the critical 
disorder reported for the Anderson model. For the simple cubic lattice nu- 
merical estimates of quantum site-percolation thresholds range from 0.4 to 
0.5 (sec [46] and references therein). In Figs. 6-8 we present data for p ty which 
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shows that p q > p c . In fact, within numerical accuracy, we found p ty = for 
p = 0.337 > p c . 

To get a more detailed picture we calculate the normalized typical DOS, 
Pty/pme, m the whole concentration-energy plane. The left panel of Fig. 10 
presents such kind of phase diagram of the quantum percolation model. This 
data also supports a finite quantum percolation threshold p q ~ 0.4 > p c (see 
also [36, 46, 47, 48]), but as the discussion above indicated, for E = and 
E = ±t the critical value p q {E) = 1, and the same may hold for the set of 
other 'special' energies. The transition line between localized and extended 
states, p q {E), might thus be a rather irregular (fractal?) function. On the 
basis of our LDOS distribution approach, however, we are not in the position 
to answer this question with full rigor. 

Finally let us come back to the characterization of extended and local- 
ized states in terms of distribution functions. The right panel of Fig. 10 dis- 
plays the distribution function, F[pi(E)], for four typical points in the energy- 
concentration plane, corresponding to localized, extended, and anomalously 
localized states, respectively. The differences in F[pi(E)] are significant. The 
slow increase of F[pi(E)] observed for localized states corresponds to an ex- 
tremely broad LDOS-distribution, with a very small most probable (or typi- 
cal) value of Pi(E). This is in agreement with the findings for the Anderson 
model. Accordingly the jump-like increase found for extended states is related 
to an extremely narrow distribution of the LDOS centered around the mean 
DOS, where p mc coincides with the most probable value. At E = and low 
p, the distribution function exhibits two steps, leading to a bimodal distri- 
bution density. Here the first (second) maximum is related to sites with a 
small (large) amplitude of the wave function - a feature that substantiates 
the checkerboard structure discussed above. For higher p, where we already 
found a reduced spectral weight of the central peak in p mc , also the two step 
shape of the distribution function is less pronounced. Therefore we may argue, 
the increase in weight of the central peak for lower p is substantially due to 
the checkerboard states. 

Having a rather complete perception of the physics in 31?, let us now come 
to the 2D case, for which the findings in literature are more controversial. 

3.2 2D Site Percolation 

Although the main characteristics of the 3D site-percolation problem, e.g., 
the fragmentation of the spectrum, persist in 2D, there exist some particu- 
larities and additional difficulties. In particular, the existence of a quantum 
percolation threshold 1 > p q > p c = 0.592746 is even less settled than in 3D. 
Published estimates range from p q — p c to p q = 1 (see [46, 49] and references 
therein). This uncertainty is due to the extremely large length scales, on which 
localization phenomena take place in 2D, a fact well-known for the standard 
Anderson model. Furthermore the special characteristics of the band center 
states seem to be of particular importance [44, 50] . 
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Local Density of States 



Figure 11 shows the mean and typical DOS of the 2D quantum percolation 
model, calculated by the LD approach. For large p, p mc clearly resembles the 
2D DOS shape for the ordered case. For these parameters pt y nearly coincides 
with p me except for the band center, where pt y shows a strong dip. If we 
reduce the occupation probability, the spikes at special energies appear again 
(see Sect. 3.1), with most spectral weight at E — 0. The weight of the central 
peak (9% close to p c ) is reduced as compared to the 3D case. 

In order to obtain reliable results for the infinite system, we examine the 
dependency of p ty on the system size, for fixed Na. Here, we find a charac- 
teristic difference between large and small p. Whereas for large p, above a 
certain system size, p ty is independent of N, it continuously decreases for low 
occupation probabilities. This behavior is evocative of extended and localized 
states, respectively. Taking a look at the underlying distribution functions, 
we find a similar situation as in the 3D case. At E = 0, the two level distri- 
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Fig. 11. Upper panels: Mean (solid line) and typical (dashed lines) DOS for the 
2D percolation model for p — 0.7 (left) and p = 0.9 (right). The dashed lines 
correspond to the typical DOS for different system sizes (from top to bottom) N = 
114 2 ,400 2 ,572 2 , 1144 2 for p = 0.7 and N = 100 2 , 350 2 , 500 2 , 1000 2 for p = 0.9. 
Again, we adapt a, such that Na = 45. The curves are ensemble averages over 
1000 (500 for the largest systems) individual LDOS spectra for different cluster 
realizations and sites. Lower panel: Distribution function of the LDOS for three 
characteristic energies: E = (solid lines), E = 0.5t (dashed-dotted lines), and 
E — 2.5t (dashed lines). To underline the finite size scaling, we show the results for 
two system sizes, 350 2 (400 2 ) marked with circles, and 500 2 (572 2 ) with squares 
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bution evolves, indicating the checkerboard structure of the state. Away from 
the special energies, the distribution function exhibits the shape characteris- 
tic for extended and localized states. This behavior exposes when we compare 
F[pi(E)] for different system sizes. Whereas for extended states the distribu- 
tion function is insensitive against this scaling, it shifts towards smaller values 
for localized ones. 

Although these results are quite encouraging, one aspect deserves further 
attention. If we try to calculate the LDOS distribution at a given energy E, 
due to the finite resolution of the KPM, it will also contain contributions from 
states in the vicinity of E. Thus taking the fragmentation of the spectrum 
into localized and extended states seriously, the LDOS distribution within 
this artificial interval will contain contributions of each class of states. 

For practical calculations, this causes no problems, as except for the most 
pronounced peaks, the probability of finding a state which is localized on one 
of the geometries like in the right panel of Fig. 8 drops very fast with its 
complexity. 

Time Evolution 

The expansion of the time evolution operator U(t) = c~ iHt in Chcbyshev 
polynomials allows for a very efficient method to calculate the dynamics of a 
quantum system. We may profit from this fact by calculating the recurrence 
probability of a particle to a given site, Pr(t), which for r — > oo may serve as 
a criterion for localization [27, 2]. While for extended states on the spanning 
cluster Pr(t — ► oo) = 1/N, which scales to zero in the thermodynamic limit, 
a localized state will have a finite value of Pr(t) as N — > oo. The advantage 
of considering the time evolution is that in general the initial state is not 
an eigenstate of the system and therefore contains contributions of the whole 
spectrum. This allows for a global search for extended states and a detection 
of p q . 

Let us briefly outline how to calculate the time evolution of the system 
by means of Chebyshev expansion [51, 52]. Of course, as a first step, the 
Hamiltonian has to be rescaled to the definition interval of the Chebyshev 
polynomials, leading to 



A I 



U(t) = e 
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The expansion coefficients c k are given by 



l 



Ck = 



f T k (x)e' [ax 

J 7T\/i — X 2 



dx = (-i) k J k (aT) , 



(17) 



-l 



where J k denotes the Bessel function of the first kind. Due to the fast asymp- 
totic decay of the Bessel functions 
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the higher-order expansion terms vanish rapidly. Thus we do not need addi- 
tional damping factors as in the Chebyshev expansion of the LDOS, but may 
truncate the series at some order without having to expect Gibbs oscillations. 
Note that the expansion order necessary to obtain precise results will surely 
depend on the time step used, normally M ~ 2at will be sufficient. Thus, for 
optimum performance of the algorithm, we have to find a suitable compromise 
between time step At and expansion order M. Anyhow, for reasonable M, 
the method permits the use of larger time steps compared to the standard 
Crank-Nicolson algorithm. 

Having this powerful tool at hand, we are now in the position to calculate 
how ip(Ti,r) evolves on Aoc in time. The 'natural' time scales of the system 
are given by the energy of the nearest-neighbor hopping element r = 1/t, 
describing one hopping event, and the time a particle needs (in a completely 
ordered system) to visit each site once, To = Nt . As initial state, we prepare 
a completely localized state, whose amplitude vanishes exactly, except for 
two neighboring sites, where it has amplitudes a and \fl — a 2 , respectively. 
For this state we can calculate the energy E = 2toVl — a 2 and, choosing a 
appropriately, we may continuously tune E e [—t, t\. Taking into account more 
complicated initial configurations of occupied sites (see right panel in Fig. 8) 
we may also adjust higher energies. For each starting position, however, the 
local structure of A x limits the possible configurations. 

In Fig. 12 we compare the time evolution of a state for high and low 
occupation probability p, for which qualitatively different behaviors emerge. 
For p — 0.65, the wave function is localized on a finite region of the cluster. 
Following the time evolution up to very long times (> 100T ) we demonstrate 
that this is not just a transitional state during the spreading process of the 
wave function, but true 'dynamical' localization. This behavior is in strong 
contrast to p — 0.90, where the state spreads over the whole cluster within 
a short fraction of this time (0.17b). For any fixed time there are some sites 
with slightly larger amplitudes (the darker dots in the last image of the time 
series). Those are due to contributions from localized states which are also 
present in the initial state. However, as the wave function extends over the 
whole cluster, for large p not all states in the spectum may be localized. 
Since the time evolution of initial states at p — 0.65 and p = 0.90 behaves in 
such a different manner, we conclude that there exists a quantum percolation 
threshold in between. 

Consequently, we have a method at hand, capable to visualize the dynam- 
ical properties of localization in quantum percolation [53] . 
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Fig. 12. Time evolution of a state at E — 0.5t on the spanning cluster of a N = 
200 2 lattice for p — 0.65 (left) and p = 0.90 (right). We show the normalized 
square of the wave function amplitude at the different sites N\ip(ri, r)\ 2 . Due to 
this normalization, for an extended state which is evenly spread over all sites of 
this quantity is equal to unity 
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4 Percolative Effects in Advanced Materials 

Current applications of quantum percolation concern e.g. transport properties 
of doped semiconductors [44] and granular metals [54] , metal-insulator tran- 
sition in 2D n-GaAs heterostructures [55], wave propagation through binary 
inhomogeneous media [56], superconductor- insulator and (integer) quantum 
Hall transitions [57, 58], or the dynamics of atomic Fermi-Bose mixtures [59]. 
Another important example is the metal-insulator transition in perovskite 
manganite films and the related colossal magnetoresistance (CMR) effect, 
which in the meantime is believed to be inherently percolative [60]. Quite 
recently (quantum) percolation models have been proposed to mimic the 
minimal conductivity in undopcd graphcne [61]. In doped graphene, in the 
split-band regime, an internal mobility edge might appear around the Fermi 
energy by introducing impurities [62]. Moreover, geometric disorder is shown 
to strongly affect electronic transport in graphene nanoribbons [63]. In the 
remaining part of this paper we exemplarily investigate two specific random 
resistor network models to describe qualitatively charge transport in graphene 
sheets and bulk manganites. 

4.1 Graphene 

Due to its remarkable electronic properties, recently a lot of activity has been 
devoted to graphene, the atomic mono/bi-layer modification of graphite [64, 
65, 66]. Especially the gapless spectrum with linear dispersion near the Fermi 
energy and the possibility of continuously varying the charge carrier density 
(from n- to p- type) by means of applying a gate voltage are of technological 
interest. In view of possible applications, like graphenc-based field effect tran- 
sistors, it is highly desirable to know how these characteristics change in the 
presence of disorder, inherent in any prepared probe. Therefore much work 
has been dedicated to study possible localization effects due to the presence 
of disorder (cf. [67] and references therein) . 

The extraordinary electronic structure of graphene results in unusual 
transport properties. In this material a finite minimal conductivity is ob- 
served, which might be attributed to a mesoscopically inhomogeneous density 
of charge carriers [68, 69, 70], caused by spatially varying charge trapping on 
the substrate. To describe the influence of these charge inhomogeneities on the 
transport properties, percolative random resistor networks (RRN) have been 
proposed [61]. Following this line, we apply the LD approach to a minimal 
model [71] that can be constructed in generalization of the 2D percolation 
model described in Sect. 3.2. 

Let us consider a 2D lattice on which the sublattices represent regions 
of different charge carrier concentrations. These regions (sites) are randomly 
connected with each other (left panel of Fig. 13). The hopping probability 
for such links (to some next nearest neighbors) is assumed to be much higher 
than for direct nearest neighbor hopping events. The later ones are reduced 
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by the leakage, A, as compared to the others. To examine the influence of 
anisotropic hopping, we model the corresponding RRN by the Hamiltonian 

H = ~4 A ^ (<4 C « + c ln C e) + ( 7 ?™ C mC„ e + (1 - Vm)4 c n) + 

(19) 

(Vmctcn + (1 - J7 m )<4c„ 1 



0(1(1 



+ H.c. , 



where e(ast) and n(orth) denote the nearest neighbors of site m in x- and 
y-direction and ne is the next-nearest neighbor in 'north-east' direction. We 
assume PBC. The random variable n m € {0, 1} determines which diagonal in 
a square is connected (right panel in Fig. 13). Shifting the expectation value 
of the {?7 m }, we may adjust the anisotropy of the system p = (r) m ). While the 
link directions are isotropically distributed for p = 0.5, increasing (decreasing) 
p generates a preferred direction of hopping, favoring stripe-like structures. 
Note that in the limit of vanishing leakage this model is equivalent to the 2D 
percolation model discussed in Sect. 3. 

First representative results for the RRN model are presented in Fig. 14. In 
particular it shows the influence of finite leakage and anisotropy. In contrast 
to the 2D percolation model, where the DOS spectra are completely symmet- 
ric, the inclusion of next-nearest neighbor hopping causes an asymmetry that 
grows with increasing A. For large p, the mean DOS is evocative of the ID 
DOS, except for the multitude of spikes, which we can attribute again to local- 
ized states on isolated islands. In the isotropic, low-leakage case, a vanishing 
pt y suggests that all states are localized. Either increasing p or A leads to a 
finite value of pt y . But even at p = 0.90 this effect is marginal for small A, thus 
presumably no extended states exist. Increasing the leakage results in a finite 
Pty for E > also at p — 0.5. This feature becomes even more pronounced 




Fig. 13. Left panel: One particular cluster realization of the RRN model at p = 0.5. 
Right panel: Visualization of the generation rule for the RRN model 
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at high anisotropy, which indicates the existence of extended states for these 
parameters. 




Fig. 14. Mean (solid line) and typical (dashed line) DOS for the RRN model. 
Results are shown for N = 500 2 and Na = 45 for different anisotropy p and leakage 
A. For illustration of the finite size dependence, pt y for a N = 350 2 sites system is 
also included (dotted lines). The results base on the calculation of fOOO individual 
LDOS spectra for different cluster realizations and sites 



In Fig. 15 we present some characteristic eigenstates of this model for a 
fixed realization of disorder. These results support the conclusions drawn from 
the typical DOS. For the isotropic, low-leakage case we find a clearly localized 
state with internal checkerboard structure. Increasing the anisotropy, many 
states are extended in one direction (cf. the lZ?-like shape of p m c), while 
localized in the other one. The leakage has a more drastic effect on the nature 
of the states, as they are extended in some sense for both values of p. The 
amplitudes on different lattice sites fluctuate over several orders of magnitude, 
however, explaining the reduced value of pt y as compared to p mc (cf. Fig 14). 

Due to the simplicity of the model, these results are surely not suitable 
to be compared to real experimental transport data, but can be seen as a 
first step towards an at least qualitative understanding of the extraordinary 
transport properties in graphenc. In any case, also here the LD approach may 
serve as a reliable tool to discuss localization effects. 
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p=0.50 ; M).01 ; E=0.000 p=0.50 ; M).1 ; E=1.402 




Fig. 15. Characteristic eigenstates of the RRN model on a N = 100 lattice at 
different anisotropies p, leakage rates A, and energies E. Shown is the normalized 
occupation probability iV|(i|n)| 2 



4.2 Doped CMR Manganites 

The transition from a metallic ferromagnetic (FM) low-temperature phase to 
an insulating paramagnetic high-temperature phase observed in some hole- 
doped manganese oxides (such as the perovskite family Lai_ x [Sr, Ca] x Mn0.3) 
is associated with an unusual dramatic change in their electronic and mag- 
netic properties. This includes a spectacularly large negative magnetoresistive 
response to an applied magnetic field (see Fig. 16, left panel), which might 
have important technological applications [72]. 

Recent experiments indicate the coexistence of localized and itinerant 
charge carriers close to the metal-insulator transition in the FM phase of CMR 
manganites. Above T c the activated behaviour of the conductivity [73] as well 
as the structure of the pair distribution function (PDF) [74] indicate the for- 
mation of small polarons, i.e., of almost localized carriers within a surrounding 
lattice distortion. Interestingly these polarons continue to exist in the metallic 
phase below T c , merely their volume fraction is noticeable reduced. For the 
coexistence of conducting and insulating regions within the metallic phase 
different scenarios were discussed, which relate the metal-insulator transition 
to phase separation [75] and percolative phenomena [76, 77]. In particular 
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microscopic imaging techniques, like scanning tunneling spectroscopy [78, 79] 
or dark- field imaging [80], seem to support the latter idea. 

In previous work [81] we addressed this problem theoretically. We proposed 
a phenomenological mixed-phase description which is based on the competi- 
tion of a polaronic insulating phase and a metallic, double-exchange driven 
ferromagnetic phase, whose volume fractions and carrier concentrations are 
determined self-consistently by requiring equal pressure and chemical poten- 
tial. In more detail, we assume that the resistivity of the metallic component 
is proportional to the expression 

derived by Kubo and Ohata [82], which associates p with the fluctuation of 
the double-exchange matrix element caused by the thermal spin disorder. Here 
5 = 3/2 is the localized spin formed by the t 2g electrons of the manganese. 
Both, 

SB s [z] /,„„ „, , (S+l)z , z\ 5 + 1 
9sW = (25TIF ( <2S + 2) C °' h - C ° th 2sj + 2STT (21) 

wW-lf^ + jg^TOU.^.jB.W, (22) 

exhibit a magnetic field dependence, where B$ [z] — [(25+ 1) coth ^ 2S ^ Z — 
coth . The resistivity of the insulating component is assumed to match the 
resistivity of the high-temperature phase, which in experiment is well fit by 
the activated hopping of small-polarons [73]. Hence, the resistivities of the 
two components are given by, 

P U) = J^y (ps[5(A + A cxt )] + Pmin ) and = P s[S\ cxt ] , 

(23) 

where e p is the polaron binding energy, (3 — 1/ksT the inverse temperature, A 
is the inner Weiss field, A cxt ist the external magnetic field and the prefactors 
A and B as well as the cut-off p m ; n are free model parameters which could 
be estimated from experimental data. Then, at a given doping level x, i.e. 
chemical potential /x, the resulting carrier concentrations and x {p) in the 
coexisting regions define the two volume fractions by the equations 

x = p(fi x m + p {p) x {p) and + = 1 , (24) 

(see [81]). The resistivity of the whole sample, which may consist of an in- 
homogeneous mixture of both components, is calculated on the basis of a 
RRN. More precisely, we choose nodes from a cubic lattice which belong to 
the metallic component with probability and to the polaronic component 
with probability p( p \ Each of these nodes, which represent macroscopic re- 
gions of the sample, is connected to its neighbours with resistors of magnitude 
pit) or p(p) ^ respectively. 
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Fig. 16. Experimental (left) and theoretical (right) results for the temperature- 
and magnetic field-dependence of the resistivity in doped CMR manganites. Exper- 
imental data are taken from Ref. [83]. For a more detailed discussion see also [81] 



Inserting the volume fractions and carrier concentrations from the mixed- 
phase model we obtain the resistivities shown in Fig. 16 (right panel). The 
jump-like behaviour of the resistivity originates to a large degree from the 
changing volume fraction of the metallic component, which can cross the per- 
colation threshold. However, the conductivity of the component itself as well 
as its carrier concentration strongly affect p for T <T C . An external magnetic 
field causes a reasonable suppression of p, i.e., a noticeable negative magne- 
toresistance. Compared to the real compounds the calculated effect is a bit 
weaker. Nevertheless, in view of the rather simple model for the conductivity 
the agreement is quite satisfactory. More involved assumptions, e.g. an affinity 
to the formation of larger regions of the same type in the sense of correlated 
percolation [84] would naturally affect the resistivity of the system and its 
response to an external field. 



5 Conclusions 

In this tutorial we demonstrated the capability of the local distribution ap- 
proach to the problem of quantum percolation. In disordered systems the 
local density of states (LDOS) emerges as a stochastic, random quantity. It 
makes sense to take this stochastic character seriously and to incorporate the 
distribution of the LDOS in a description of disorder. Employing the Ker- 
nel Polynomial Method we can resolve with very moderate computational 
costs the rich structures in the density of states originating from the irregular 
boundary of the spanning cluster. 
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As for the standard Anderson localization and binary alloy problems the 
geometrically averaged (typical) density of states characterizes the LDOS dis- 
tribution and may serve as a kind of order parameter differentiating between 
extended and localized states. For both 2D and 3D quantum site percola- 
tion, our numerical data corroborate previous results in favor of a quantum 
percolation threshold p q > p c and a fragmentation of the electronic spectrum 
into extended and localized states. At the band center, so called anomalous 
localization is observed, which manifests itself in a checkerboard-like struc- 
ture of the wave function. Most notably, monitoring the spatial evolution of 
a wave packet in time for the 2D case, we find direct evidence for 'dynami- 
cal' localization of an incident quantum particle at p = 0.65 > p c , while its 
wave function is spread over the percolated cluster for p = 0.9. This finding 
additionally supports the existence of a quantum percolation threshold. 

Without a doubt quantum percolation plays an important role in the trans- 
port of several contemporary materials, such as 2D graphene or 3D manganese 
oxides. To close the gap between the study of simple percolation models and 
a realistic treatment of percolative transport in these rather complicated ma- 
terials will certainly be a challenge of future research. 

We thank A. Alvermann, F. X. Bronold, J. W. Kantelhardt, S. A. Trug- 
man, A. Wcifie, and G. Wcllein for valuable discussions. This work was sup- 
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work for Technical/Scientific High-Performance Computing in Bavaria (KON- 
WIHR), and the Leibniz Computing Center (LRZ) Munich. 
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